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§ ■ Abstract 

We simulate model systems of charged spherical particles in their counterion 
i-rt I solution and measure the thermodynamic pressure and the pair distribution 

C ' function from which we derive effective potentials of mean force. For a system 

with only electrostatic and hard core interactions, we investigate the effective 
potential between two like-charged spheres in divalent counterion solution as 
a function of concentration. We find a strong attractive interaction for high 
^^ I concentration and a global repulsive effective interaction for dilute systems. 

l/^ ' The results indicate a first order phase transition in sphere-counterion density 

^^ ■ as a function of global concentration and the effective sphere-sphere potentials 

Q . in the dilute (solvated) regime suggest significant density fiuctuations due to 

^N I short range local minima in the effective energy surface. Our results arise 

from a minimal approach model of several recent experiments on polystyrene 
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C^ ' latex particles in monovalent counterion solution. 
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The effective interactions between charged macromolecules in countercharged aqueous 
media have been assumed to be governed mainly by simple screened Coulomb behavior 
leading to overall repulsion between like charged objects |]l]-0]. This notion arises primarily 
from mean-field analyses, such as Poisson-Boltzmann equations, of the statistical counte- 
rion distribution, and leads to the well known Derjaguin-Landau-Verwey-Overbeek (DLVO) 
purely repulsive effective potential between charged colloidal spherical macroions PJ^ . How- 
ever, recent experiments have indicated that the effective interaction can indeed be attractive 
in dense suspensions of, e.g., polystyrene latex particles in aqueous solutions with monova- 
lent counterions ]5|-P^. These experiments, performed with very low salt concentrations, 
have observed inhomogeneous colloid densities under certain conditions and almost perfectly 
homogeneous densities under others. Some of the experiments observe the attractive po- 
tentials only when the colloidal spheres are suspended between glass walls [^,0 and it is 
consequently a speculation whether the effective attraction is mediated by these walls or 
if the walls merely provide a confinement of phase space such that the attraction is en- 
hanced. Recent advances for investigating the effective interactions have been made |16|] by 



performing direct numerical simulations of charged colloidal spheres and their monovalent 
counterions in confined space and only repulsive interactions were found. 

As the present state of theoretical predictions is almost entirely limited to the trivial 
repulsive behavior between fractionally screened like-charged objects, we have carried out 
simulations of a simple system which exhibits strong effective short as well as long range 
attractions between like-charged spherical objects in their neutralizing counterion solution. 
This is intended as helping to construct a successful theoretical approach describing the 
complex phase behavior of charged colloidal suspensions. The purpose of this paper is 
therefore not to give a direct interpretation of the recent experiments on polystyrene latex 
particles, but instead point to a simple system exhibiting a phenomenon that a theoretical 
approach should reproduce. 

We present simulations of charged colloidal spheres and their counterions confined in 

a cubic box volume V = L^, where L is the box length in each direction, with periodic 

boundary conditions applied such that a finite bulk concentration is simulated. Finite con- 

gj centration is important since Coulomb interactions between spherical objects (~ 1/r) are 

§ inferior to entropic repulsion (~ In r) at large distances. The equilibrium distance between 

^ any charged spherical objects at finite temperature is consequently divergent in the limit 

4 of N/V -^ 0, where A^ is the number of spherical objects; thus, the effective interaction 
g is always repulsive in the dilute limit. We will focus on a system with divalent counteri- 

5 ons in order to clearly demonstrate the importance of concentration and correlation effects, 
"3 which are absent in mean field treatments leading to, e.g., the DLVO potential. We demon- 
^ strate that the effective potential between two spheres, as derived from the pair distribution 
3 function, can become attractive in finite concentration. Our simulations are performed as 

a minimal approach in that only long-range Coulomb interactions and short-range volume 
exclusions contribute to the evolution of the system. 

The model under consideration is given by the Overdamped Langevin equation of motion 
for the ith particle: 



Uiri = -ViE + ni{t) , (1) 

where fj is the normalized coordinate of the ith particle moving with a normalized friction 
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coefficient z/j (note that neither the friction coefficient nor the time normahzation has any 
importance for the equihbrium properties of the system). The total normahzed energy E of 
the system is given by, 
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E = Y.Y.[—yir^,) + WUr)] , (2) 

«=1 j>i 



where 



WUr) 



, {n, - r%f > 2A,/B,, 

Aij and Bij are the Lennard- Jones volume exclusion parameters. Notice that we only con- 
sider the short range repulsive part of the Lennard- Jones potential. A^ is the total number 
of coordinates (spheres and counterions) and eg^ is the charge of the ith particle. The 
normalized Coulomb interaction between charges eqi and eqj at a normalized distance rij 
is represented by ^^V{rij), with dielectric constant e (for details on exact summation of 
Coulomb interactions in orthorhombic systems, see Ref. fl^)- The thermal equilibrium with 
the surroundings, i.e., with the water which is not considered except as a homogeneous di- 
electric, is maintained by the noise term ni{t) which is linked to the dissipation- fluctuation 
theorem [|I^] by the uncorrelated white Gaussian distribution, 

{ni{t))=0 , {ni{t) ■ nj{t')) = 6ui^6ij6{t - t') , (4) 

where 5ij is the Kronecker delta function and 6{t — t') the Dirac delta function. Boltzmann's 

constant is denoted fc^, the temperature is T, and the energy is normalized to ii^o=l kcal/mol. 

The factor, 6, is due to n being a three dimensional vector. Throughout this paper we have 

chosen divalent counterions, sphere charges of qsp = —10, Aij ^ 21-10^, Bij ^ 30, and 

^ £ = 80 in order to mimic bulk screening of water. The hard core parameter, r^-, is chosen as 

ON r% = for counterion counterion interactions, r° =5.5 for counterion sphere interactions, 

X and r°- = 11 for sphere sphere interactions. Each simulation is initiated with random, non- 

4 overlapping positions of the counterions and spheres and a large number of time steps, of 
Oh order 10*, are discarded before accumulating the pair distribution density function, p(r), 

5 between spheres. In equilibrium, this density is directly related to the Boltzmann factor. 



> p{r)r~.exp{-W{r)/kBT) , (5) 



where W{r) is the effective potential of mean force (PMF). Inverting this relationship and 
normalizing the potential to the value at maximum distance r = L/2, we obtain the derived 
effective PMF from, 

»nr) = -*.rin(^). (6) 

Figure 1 displays the PMF between two spheres and their accompanying divalent coun- 
terions in a simulation box of volume V = L^. The friction coefficients are set to z/^ = 1 for 
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counterions and Usp = 15 for spheres with the accompanying reduced time step dt = 0.005. 
The resulting sphere-sphere PMF, derived over ~10^ time steps, reveals strongly attractive 
behavior for small L (lower curves), but as L is increased, the available volume increases, 
and the minimum characteristic for attraction lifts above zero (normalized to the potential 
at distance L/2) and eventually disappears. The minimum of the PMF, Wmin, is shown 
in Fig. 2 as a function of the simulated volume. Here we observe how the minimum (rel- 
ative to the value at L/2) of the effective potential depends on the volume. At a volume 
«i (120A)^ «i 2 ■ lO^nm^ the potential at L/2 becomes the global minimum, which clearly in- 
dicates a transition from close packed spheres being preferred in small volumes to separated 
spheres being statistically preferred at larger volumes. 

However, as is obvious from Fig. 1, the spheres can be attracted to each other at short 
distances even if the global PMF minimum is at L/2. This feature arises from the much 
stronger interaction that exists between a counterion and two closely spaced spheres when 
compared to the interaction between a counterion and a single sphere. This then leads 
to a compact state of the two spheres with their counterions spatially correlated around 
them. The barrier height (PMF) in energy going from the compact state to the maximum 
separation {L/2) is shown in Fig. 3, where we find a significant barrier, i.e., a local minimum 
for the compact state, well into the volume regime where the global minimum of the PMF 
is at maximum separation. 

The PMF potentials shown in figure 1 are in strong contrast to the DLVO potential, 
which is given by MjM 






W^DLVo(r) = -^qL ^ e--'-^ + W,,{r) , (7) 



where a is the radius of the spheres and the hard core interaction between two spheres is 
given by Whc- The two characteristic length scales, A^ and n, are given by, 

(8) 
(9) 
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where M is the number of counterions and qd is the counterion valency. The DLVO potential 

is obviously repulsive and monotonically decaying with distance. 

Studying the PMF gives information about the spatial correlations within the sample 
^ simulation, but it does not reveal the overall nature of global attraction versus repulsion of 

> an ensemble. In order to study if the system tends to collapse or expand, we have calculated 

^ the thermodynamic pressure. 
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where JF is Helmholtz's free energy and (■ ■ ■) denotes a temporal average. The first term on 
the right hand side is the thermal contribution to the pressure, the second is the contribution 
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from the virial of all hard core interactions (^, and the last is the electrostatic contribution, 
Uei being the total electrostatic energy of the system. 

Figure 4a shows the total pressure, px, versus the volume for one simulated sphere and 
its countreions in a cubic box. In agreement with the PMFs visualized in figure 1, we find 
an effective attraction (negative pressure) for high concentration of spheres and a repulsive 
effective interaction (positive pressure) for dilute systems. This behavior of pressure as a 
function of concentration indicates a first order phase transition from a solvated state of 
spheres at low concentration to a collapsed state of spheres at high concentration. However, 
the global pressure cannot reveal the short range local minima in the PMFs observed in 
figure 1. These local minima, which persist far into the dilute limit, will result in local 
density fluctuations of the solvated spheres. Note that the negative pressure is solely due 
to the electrostatic interactions since both the kinetic and the hard core contributions to 
the total pressure are always positive (see figure 4b, where the three individual pressure 
contributions are shown). The total pressure will, in the dilute limit, be entirely dictated by 
the thermal contribution (ideal gas) since the finite temperature will separate all spherical 
particles and thereby make the hard core and the electrostatic contributions negligible. 

Also the pressure is in contrast to the DLVO analysis. Since the DLVO potential is purely 
repulsive, one can only expect positive pressures from such a treatment. The reason for the 
possibility of negative pressures in the "all particle" model lies in the overall charge neutrality 
of the system. It is clear that a charge neutral system composed of oppositely charged 
components is self- contracting at low temperatures (e.g., the sodium- chloride crystal) due 
to strong spatial correlations between the charged objects. However, the DLVO analysis 
only accounts for the counterion charge if it contributes to the screening of a sphere, and 
thus, the DLVO treatment does not maintain the overall charge neutrality of the system. 

We have simulated a system of charged spheres in divalent counterion solution and 
demonstrated a very strong attractive effective potential between the like-charged spheres. 
The origin of the attraction lies in the strong correlation effects between counterions when 
the spheres are in a compact conformation. This simple result is in contrast to the effective 
(and always repulsive) DLVO potential, which is developed without considering the coun- 
terions as particles. However, our results are in agreement with other recent simulations of 



ON like-charged rods in divalent counterion solution |]I9| and theoretical considerations regard- 

X ing like-charged plates |^. Our simple simulations demonstrate that a successful theory 

4 explaining the complicated phase diagram of complex charged fluids needs to account for the 
g correlated behavior of counterions, the global charge neutrality leading to the tensile pres- 

5 sure in figure 4, and the relationship between the sphere radius and the simulated volume 
"3 responsible for the statistical distribution of counterions. 

^ It is important to emphasize that our simulations are not meant to characterize real 

3 experiments quantitatively. The minimal approach model presented in this paper cannot 

account for any atomic or molecular detail that may be important for experiments. For 
example, including the effects of the solvent (water) as friction, noise, and bulk screening is 
a dramatic simplification of the true discrete nature of water molecules, which exhibit their 
discrete nature within the hydration shell. Also, neglecting the complexity of the dielectric 
mixture of water and spheres prevents our model from correctly simulating the electrostatic 
problem. Our simulations therefore address the validity of present theories, such as DLVO, 
which attempt to describe the minimal approach model described in this paper. 
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FIG. 1. Potentials of Mean Force for different simulated volumes. The potentials are set to 
zero at maximum separation, L/2, and the different potentials are off-set vertically in order to 
better distinguish the cases. Potentials of Mean Force are shown for (from the top): L = 200A, 
L = 150A, L = 120A, L = lOOA, L = 80A, L = 70A, L = 60A, and L = 50A 
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FIG. 2. Minimum of the potentials of mean force shown in Fig. 1 versus the simulated volume, 
V = L^. 
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FIG. 3. Energy barrier height from the compact state (r ~ 16 in figure 1) to the local PMF 
maximum (r ~ 22 in figure 1) as a function of the simulated volume. 
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FIG. 4. (a) Magnitude of thermodynamic pressure of a system of one sphere and its counterions 
in a cubic box versus simulated volume. The inset shows the same data in a log-log representation, 
(b) Magnitude of the three contributions to the pressure shown in (a); thermal (dotted); hard core 
(sohd); and the negative electrostatic (dash-dot). 
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